function cc=seeiwsv(priorV,priorT,ndr); 
% ===========================================
% SEEIWSV.m
% function cc=seeiwsv(priorV,priorT,ndr); 
% Plot and provide percentiles of an Independent IW 
% distribution as parametrized for the SV projects, i.e. 
% degrees of freedom = priorT 
% scale matrix       = priorV*priorT 
% 
% STD is finite only for priorT > 4 
% NDR is optional and default is 50,000
% 
% cc is the cell printed onscreen 
% 
% Note: this is only for uncorrelated IW 
%
% Alejandro Justiniano April 25 2007 
% =============================================
nx=size(priorV,1); 
[junk,check]=chol(priorV*priorT); 
if check~=0; error('PriorV*PriorT must be a PD matrix'); end; 
if nargin < 3 | isempty(ndr); 
    ndr=50000; 
end 
vdraws=zeros(nx,ndr); 
vinv=eye(nx)/priorV*priorT; 
for ii=1:ndr; 
    shock=mvnrnd(zeros(nx,1),vinv,priorT);
    %vdraws(:,ii)=1./diag((shock'*shock));
    vdraws(:,ii)=diag(eye(nx)/(shock'*shock)); 
end 
vdraws=vdraws'; 
% Theoretical and Simulated STD 
mean_mat=zeros(nx,2); 
std_mat=zeros(nx,2); 
med_mat=zeros(nx); 
pvec=round(ndr*[0.01 0.05 0.1 0.5 0.9 0.95 0.99]); 
perc_mat=zeros(nx,7); 
vdraws=sort(vdraws,1); 
nrow=round(nx/2); 
namecell=fnumcell('V',[1:nx]); 
figure; 
ii=1; 
for ii=1:nx; 
    subplot(nrow,2,ii); 
    hist(vdraws(pvec(2):pvec(end-1),ii),100); 
    title([namecell{ii},' priorV:',num2str(priorV(ii,ii))]); 
    mean_mat(ii,1)=priorV(ii,ii)*priorT/(priorT-2); 
    mean_mat(ii,2)=mean(vdraws(:,ii));
    if priorT > 4
    std_mat(ii,1)=sqrt( 2*(mean_mat(ii,1)^2)/(priorT-4) ); 
    else 
        std_mat(ii,1)=inf; 
    end 
    std_mat(ii,2)=std(vdraws(:,ii)); 
    med_mat(ii)=median(vdraws(:,ii)); 
    perc_mat(ii,:)=vdraws(pvec,ii); 
    vline(med_mat(ii)); 
    axis tight; 
end 
suptitle_withpatch(['Hist for 5-95 & priorT:',num2str(priorT)]); 
% Create cell with moments 
c1=crtcell([diag(priorV) mean_mat std_mat],namecell(:),...
    {'Prior V','Th Mean','Sim Mean','Th Std','Sim Std'},['Prior T:',num2str(priorT)]); 
c2=crtcell([diag(priorV) perc_mat],namecell(:),...
    {'Prior V','1%','5%','10','Med','90','95','99'},['Prior T:',num2str(priorT)]); 
cc=merge_cells(c1,c2,1); 
printcell(cc); 